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ABSTRACT 

We present Monte Carlo simulations of dusty spiral galaxies, modelled 
as bulge + disk systems, aimed to study their extinction and polarization 
properties. The extinction parameters (absorption and scattering) of dust grains 
are calculated from Mie's theory for a full distribution of sizes and materials; 
the radiation transfer is carried on for the four Stokes parameters. Photometric 
and polarimetric maps of galaxies of different optical depths, inclinations and 
bulge-to-total ratios have been produced in the B and I bandpasses. As expected, 
the effect of scattering is to reduce substantially the extinction for a given 
optical depth, in particular for what concerns the obscuration of bright bulge 
cores. For the same reason, scattering reduces also the reddening, as evaluated 
from B-I maps. On the other hand the bluing directly due to forward scattering 
is hardly appreciable. Radial color gradients are often found. A comparison with 
"sandwich" models shows that they fail dramatically to reproduce the extinction 
- optical depth relation. The degree of linear polarization produced by scattering 
is usually of the order of a few percent; it increases with optical depth, and with 
inclination (i < 80°). The polarization pattern is always perpendicular to the 
major axis, unless the dust distribution is drastically modified. There is little 
local correlation between extinction and polarization degree and there is a trend 
of increasing polarization from the B to the I band. We discuss implications and 
relevance of the results for studies of the structure and morphology of spiral 
galaxies and of their interstellar medium. 

Subject headings: dust, extinct ion — galaxies: spiral — methods: numerical- 
polarization — radiative transfer — scattering 
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1. Introduction 

For a long time, since Holmberg's first analysis (Holmberg 1958), disk galaxies, 
including the gas-rich late morphological types, have been taken as largely transparent 
systems. Although with various refinements, this point of view has been shared by the 
authors of catalogs of galaxies, such as the relatively recent Revised Shapley-Ames Catalog 
of Bright Galaxies (Sandage & Tamman 1981), and thereby adopted by the majority of 
the astronomical community. Holmberg inferred average values of the internal extinction 
on statistical grounds, that is by interpreting the correlation between inclination and mean 
surface brightness in a rather limited sample of 119 spirals of various morphologies. On the 
other hand, refined studies of the dust content and distribution of individual, representative 
galaxies are scarce. In the edge-on Sc galaxy NGC 891, Kylafis & Bahcall (1987) have been 
able to determine an overall transparency of the disk, when seen face-on, in a bandpass 
around 6000A. 

Recently, however, some studies have shaken the widespread belief in transparent 
disks (Disney, Davies & Phillipps 1989; Valentijn 1990) and the issue has rapidly become 
one of the most contended in the literature on normal galaxies, sometimes with opposite 
conclusions drawn by the analysis of the same, or similar, set of data (Valentijn 1994). 
The idea of opaque galaxies has gained in the last years a certain acceptance and has been 
shared by the recent Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 
1991). 

This debate on the optical thickness of galactic disks is relevant to a number of issues 
important for extragalactic astronomy. We mention: (i) the morphology and internal 
structure of galactic components, such as bulges, disks, lenses, bars, etc; (ii) the dark 
matter content of galaxies and clusters, and (iii) the luminosity corrections affecting 
redshift-independent measures of distances. The last point, in particular, touches on the 
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field of the large scale structure of the universe, and of the deviations from the Hubble flow. 

For these reasons, observations and data analyses are becoming more refined with 
particular reference to selection biases (Burstein, Haynes & Faber 1991), to redshift 
information (Choloniewski 1991; Giovanelli et al. 1994) and to multiwavelength imaging 
(Peletier et al. 1994; Evans 1994). At the same time the modelling has also become 
increasingly sophisticated. In fact, as it has been argued by most of the already cited 
authors, a fundamental problem in statistical evaluations of the extinction lies in the 
absence of accurate and sensitive observable indicators. Surface brightness, diameters and 
colors have often been interpreted in the framework of models so naive and unrealistic to 
yield too simplistic and sometimes wrong results. 

A realistic modelling of dusty disks is thus strongly required. Besides, if we postulate 
the possibility of large optical thicknesses and sizeable grain albedos, as it is the case for 
wavelengths shorter than 3/xm, then the effect of multiple scattering in composite (bulge + 
disk) systems, with dust interspersed with stars, must be included (e.g. Block et al. 1994). 
A remarkable progress on this track is achieved in the recent paper by Byun, Freeman & 
Kylafis (1994, hereafter BFK) who have produced simulations of disk galaxies of various 
morphology and optical thickness. Their treatment of the radiative transfer includes exactly 
the contribution of the first scattering, and approximately the higher scattering orders. 
They adopt a Henyey-Greenstein phase function and a single dust grain with average 
optical properties according to Draine & Lee (1984). 

Scattering events are also expected to introduce polarization in the otherwise 
unpolarized emission from normal stars. Polarimetric studies of normal spirals, both 
observational and theoretical ones, are scarce. In view of the current interest in the 
properties of dusty disks such scarcity is remarkable since most of the information about the 
physical properties of dust grains is locked into their polarization effects. Among the few 
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observations of this kind, we mention the detailed studies of M31, M82, M104, and of the 
Magellanic clouds; they are summarized and referenced by Wielebinski & Krause (1993). 

One of the fields where the analysis and simulation of polarization maps has been 
most effective in recent years is star formation, and, in particular, the study of protostellar 
envelopes (see Fischer, Henning & Yorke, 1994 and references therein). Here the more 
recent and detailed models are based on Monte Carlo simulations which allow for arbitrary 
distributions of polydispersive mixtures of spherical dust particles. In all these models any 
observed polarization is imputed to transfer effects on light which is otherwise emitted 
unpolarized by stars. In particular, since dust grains are usually assumed to be spherical, 
homogeneous and optically isotropic, only the polarizing effect of scattering is taken into 
account. It is widely known, anyway, that in our Galaxy the linear polarization introduced 
by the light propagation through the interstellar medium is best described as the effect 
of dichroic absorption by aligned, elongated grains (Spitzer 1978; Gledhill & Scarrot 
1989). Sometimes, this polarization mechanism has been questioned but, indeed, only for 
particular situations, where the presence of a suitable toroidal magnetic field, that is the 
main ingredient for the alignment, is dubious or improbable (Menard 1989). 

We present here the results of Monte Carlo simulations of realistic model galaxies within 
the frame of Mie's scattering theory, that is for spherical dust grains with homogeneus and 
isotropic optical properties. In our study the dust grain is not a single, average one but 
instead we allow for a mixture of different materials and a continuous distribution of size. 
The phase function is consistently computed from Mie's theory for each grain. The radiation 
transfer is exact to any scattering order and complete for the four Stokes parameters, so 
that the output can be analyzed in the form of polarimetric, as well as photometric, maps. 
Stars and dust three-dimensional distributions are the ones commonly accepted, and similar 
to those adopted by BFK. The present paper deals only with scattering events on spherical 
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grains and does not include any dichroism in the absorption itself. The polarization 
properties may, therefore, be different from those produced by needles rotating, as widely 
believed, in a galactic meridian plane. 

Although we deal in this study with the simplest of the realistic cases, that is, smooth 
distributions and spherical grains, the Monte Carlo method is ideal to probe more hostile 
situations. The effect of clumping in the obscuring dust, of non-spherical dust grains and of 
extinction by gas at short wavelengths are deferred to future papers. 



In the following we describe the model distribution here adopted to simulate the spatial 
distribution of the photon emitters (stars) and of the particles contributing to extinction 



We describe the stellar spatial distribution in a spiral galaxy by two components: a 
spheroidal bulge and a 3-dimensional disk. We neglect dark massive halos, since we are 
only interested in luminous components, and also possible small-scale inhomogeneities such 
as spiral arms and star-forming regions. 



2.1.1. Bulge 

The surface brightness profiles of spheroidal systems is usually well described by the 
i? 1 / 4 law (de Vaucouleurs 1959): 



2. A dusty galaxy model 



(dust). 



2.1. 



Stellar distributions 



I(R) = h 10 



■3.33[(,R/R e ) 1/4 -l] 



(1) 
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where R e is the equivalent radius, and I e is the surface brightness of the corresponding 
isophote; both R and R e are measured on the sky plane. For numerical purpose, it is 
convenient to have an analytical expression for the bulge luminosity density p^ b \ which can 
be obtained deprojecting the profiles. Since there is no simple, analytical expression for p^ 
corresponding to equation (|l|) (see Young 1976) we choose to adopt the Jaffe distribution 

which corresponds to the luminosity density (Jaffe 1983): 

0) 

pW(f) = ^ (2) 

(r/r b ) 2 [1 + (r/r b )] 

where f is the physical distance from the galactic center, and r b is the scale radius for wich 
we adopt r b = 1.16 R e (see Appendix [A.l.lp . 



In principle, there is no need for an external truncation radius for the Jaffe distribution 
since the volume integral converges for f — > oo; for numerical purposes, though, we have 
introduced a truncation radius rj^ ax = nr b . In fact, in absence of an external truncation, at 
large radii the bulge's light will dominate once again that of the (exponential) disk. In the 
simulations presented here, we have adopted the values R e = 1.60 kpc, n = 8. The value for 
R e is close to the one inferred for our Galaxy; the value for n is such that r$ ax corresponds 
to the point where the ratio between the bulge luminosity density and the density of the 
adopted disk is maximum (Sect. |2.1.2| ). 

2.1.2. Disk 

The assumed luminosity density distribution for the disk is (Freeman 1970) 

p«(r, z) = p { d) exp(-r/r d )Z(z/z d ), (3) 



where r and z are the cylindrical radius and the height above the equatorial plane of the 
galaxy, respectively, and the costants r d , z d are the relative scale lengths; The function 
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Z(z/zd) describes the vertical behavior of the distribution, which is assumed to be either 
exponential or of the type "sech 2 " (van der Kruit & Searle 1981; Wainscoat, Freeman 
& Hyland 1989); the results presented in this paper are for the exponential vertical 
distribution that is 

Z(z/z d ) = exp(- | z | /z d ). (4) 

In analogy with the bulge case, we have introduced a horizontal, r$ ax = mra, and a vertical 
z mlx = k z d external truncation. Our models are for r d = 4.0 kpc, z d = 0.35 kpc, which 
are similar to what is observed for the old disk population of the Galaxy; for the external 
truncations we adopt m = k = 6. 

2.1.3. bulge/total luminosity ratios 

To simulate different Hubble types, we have run models with different bulge/total 
luminosity ratios ranging from (resembling an Sd galaxy, BT00) to 0.5 (Sa, BT05). All 
our simulations are monochromatic and we have run models at the effective wavelengths 
of the B and I bandpasses, 4400 and 9000A respectively. Possible differences in the stellar 
population of the two components are taken into account by adopting a B-I color of the 
bulge larger than that of the disk. More precisely the luminosity of the bulge has been 
increased by one magnitude in the I band while the luminosity of the disk has been kept 
constant in both bands; therefore we suppose an intrinsic B-I=0 for the disk and B-I=l for 
the bulge. 

All of these prescriptions conform to those adopted by BFK. 



- 9- 



2.2. Dust properties 

The only absorbing particles included in our code are dust grains; we neglect the 
absorption above the Lyman limit due to diffuse hydrogen gas, the possible presence of 
molecules, etc. . 

2.2.1. Dust disk 

The dust is supposed to be smoothly [] distributed in a disk with the same law (eq. ||) 
as for the stellar disk 

pWfa z) = p { g) exp(-r/r g )Z(z/z g ). (5) 

where r g and z g are the horizontal and vertical dust scale lengths. In this paper we are 
mainly interested in determining the effects of a standard dust layer on the structural 
parameters and appearance of a galaxy, we will devote a following paper to study the 
detailed structure and distribution of dust as inferred from the observed light distribution in 
spiral galaxies. The Galactic parameters for the dust distribution are by far more uncertain 
than for the stars. For the computations presented here, we have assumed r g = 4.0 kpc, 
z g = 0.14 kpc and the same truncation factor as for the stellar disk. That is, while the 
radial extents of the stellar and dust disk are identical, we will always deal with dust disks 
which are vertically embedded into the star distribution. This is an important point, since 
several of the results in the following will be determined by the presence of sizeable stellar 
emission outside the obscuring layer. 

We assume that the dust physical properties (i.e. the grain size distribution and 
materials, see below) are independent of the position in the galaxy and therefore the 

1 A study of the effects introduced by dust clumping is in progress. 



absorption coefficient at wavelength A can be written as k\ = (C ext (X))p^(r, z), where 
(C ext (\)) is the extinction cross section per unit mass, averaged over the grain distribution 
and materials (see eq. f7j). 

It is common to normalize the total content of dust in the galaxy in terms of the optical 
depth along the simmetry axis of the galaxy, 



both for the exponential and sech 2 distributions. As for the stellar disk, the results here 
presented refer to vertically exponential dust disks. 



In our model the grains are supposed to be spherical and to have a size distribution 
given by the MRN model (Mathis, Rumpl & Nordsiek 1977), n(a) oc a -3 ' 5 , where a is 
the grain radius. We consider three materials: astronomical silicates, || graphite, and _L 
graphite. The numerical silicates/graphite ratio is 1:1.12, with 1/3 of the graphite having 
optical properties measured parallel (||), and 2/3 perpendicular (_L) to the c-axis. The lower 
and upper limits of the distribution are a_ = 0.005 /im and a + = 0.25 /im, irrespectively 
of the material. Very small grains and PAH have not been included due to the large 
uncertainties both in their size distribution and optical constants. 

The dielectric constants adopted are the ones given by Draine & Lee (1984), recently 
extended in the far UV and X-rays by Martin & Rouleau (1991). All the relevant 
optical properties (absorption and scattering cross section, albedo, phase function, and 
Mueller matrix) have been calculated using Mie formulae. The resulting extinction curve is 
illustrated in Fig. [[]. 




2.2.2. Dust components and size distribution 
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EDITOR: PLACE FIGURE |T| HERE. 

In our models the central optical depth of the disk is defined by Ty(0) the value at 
5500A; models at different wavelengths are still defined by Ty(0) but all computations of 
optical thickness are scaled according to our extinction curve. 



2.2.3. Comparison with Henyey-Greenstein approximation 
The average value of a quantity q, over the dust distribution, is here defined as 

= ZiW i f:+qCr t n i (a)da 

where Cf xt is the extinction cross section, the index i referes to the material, and Wi is its 
weight in the distribution; both q and C? xt depend on A 

Previous works dealing with scattering problems have often approximated the phase 
function $ with the Henyey-Greenstein (HG) formula (Henyey & Greenstein 1941), 

n 1 l -9 2 



47r(l+£ 2 -2#cos#) 3 / 2 ' 



where g is the asymmetry parameter (that is the average value of cos 6 weighted by the 
phase function); 9 and are the polar and azimuthal scattering angles, respectively. 

In Fig. ^| we show a comparison between the average phase function for the distribution 
calculated with the Monte Carlo procedure implemented in our radiative transfer code and 
the HG phase function for which a value of g averaged on the distribution according to 
equation (|?|) has been used. The differences between the two phase functions are striking (as 
an example, the B and I bands are shown): the HG approximation always underestimates 
the forward scattering by as much as 25 %, the difference being larger in the B band. 
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EDITOR: PLACE FIGURE | HERE. 

For the albedo and asymmetry parameter we obtain (uj) = 0.53 and (g) = 0.38 in the 
B-band, and (u) = 0.50 and (g) = 0.21 in the I-band. 



3. Monte Carlo simulation 

Monte Carlo methods allow the life of each photon to be followed through scattering 
and absorption processes. Here we provide a brief description of our computational scheme, 
deferring the interested reader to the Appendix for a more detailed presentation. 



3.1. Emission cycle 

We determine first the position of the photon emission point, according to the bulge 
and disk luminosities as described in Sect. |2.1| . Photons are emitted isotropically and 
unpolarized. The Stokes parameters are defined for each photon on the plane perpendicular 
to its propagation vector, following the conventions in Shurcliff (1962); the reference 
direction on such plane is the projection of the galactic z axis (i.e. the simmetry axis). 
All photons are emitted with the same intensity, or strength, that is with the first Stokes 
parameter 1 = 1. The emission cycle requires five random numbers: three for the position 



and two for the direction (see Appendix |A.1|). 



3.2. Scattering cycle 

Then we determine the position of the scattering point, i.e., the position at which the 
photon collides with a dust grain; grains are assumed to be distributed as described in Sect. 
[2.2| . This implies to determine first an optical depth r for the scattering point, and then the 
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geometrical position along the line of propagation where this r is actually reached. It turns 
out that this second step is the most consuming, in terms of computational time, of the 



whole process (Appendix |A.2|) . Once the position of the scattering event is known, we pass 
to determine the characteristics of the dust grain, according to the size and composition 
distributions (Appendix |A.3|) . Note that in principle such a sequence renders rather 
straightforward to allow for situations where the physical properties of dust vary within the 
galaxy, e.g. with galactocentric distance. These steps require two more random numbers, 
one to locate the diffusion point, the other to fix the size and composition of the scatterer. 
Knowing the radius and refractive index of the grain (together with A), it is now possible 
to compute the (polar) phase function and the albedo, and, by extracting a new random 
number, to choose the (polar) angle of scattering 9. Together with the geometro-optical 
properties of the dust grain, and within the frame of Mie's theory, 9 determines completely 
the elements of the scattering (Mueller) matrix. The matrix elements also enter into the 
expression for the azimuthal phase function, so that, by extracting a new, final random 
number we define completely direction and Stokes vector of the scattered photon (see 
Appendix Kl). 



The updated Stokes parameters are calculated via the following transformation 

( V \ 



Q' 
U' 



oc R'MR 



/ l \ 

Q 
u 



(9) 



The rotation matrix R transforms the initial Stokes parameters into a new set whose 
reference direction for linear polarization lies in the scattering plane. M is the relevant 
Mueller matrix with reference to the scattering plane, which depends only on 9 and on the 
dust properties; the matrix R' provides the final transformation of the Stokes parameters 
to the final reference direction, that is the projection of the z axis onto the plane normal 
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to the new direction of propagation. In our formulation we make use of normalized Mueller 
matrices, the absolute values of the Stokes parameter are adjusted according to the incoming 
intensity I and to the albedo of the scatterer. In summary a single scattering cycle needs 
the extraction of four random numbers: one to locate the diffusion point, one to choose the 
size and type of grain, and two for the direction of the scattered photon. The scattering 
cycle is iterated until the photon escapes the dust layer, or dies because its strength I goes 
below a threshold value (Iu m = 10~ 4 ). 



3.3. Output and performance 

The output of the code consists of Nb maps, for each Stokes parameter, where Nb is 
the number of bands of different galaxy inclination to the line of sight (Appendix |A.5|) . The 



models presented here are obtained with Nb =15, and each map has an extent of 201x201 
pixels, with a spatial resolution of 0.2 kpc; the resolution is therefore moderate both in 
space and in inclination. 

A full model (bulge+disk) with enough signal to noise ratio to distinguish a clear 
polarization pattern in the outer parts of the galaxy (with linear polarization degree of 
~ 1%) required a number of photon launches of the order of 10 8 . Anyway, photometric 
accuracies comparable to those of astronomical images are usually attainable with less 
than 10 7 launches. The code has been tested on a variety of more simple and predictable 
situations: point sources within spherical scattering envelopes, Rayleigh scattering, optically 
thin regimes, absorption only, etc.. 
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4. Results 

In summary, the models presented here were all obtained for: bulge effective radius 
R e = 1.6 kpc, truncated at 15 kpc; disk radial scale length = 4.0 kpc, truncated at 
24 kpc; disk exponential vertical scale length Zd = 0.35 kpc, truncated at 2.1 kpc; dust 
radial scale length r g = 4.0 kpc, truncated at 24 kpc; dust exponential vertical scale length 
z g = 0.14 kpc, truncated at 0.84 kpc. All of the results obtained do not depend on the 
adopted truncation factor as far as they exceed 5 scale lengths. 

We ran simulations for 6 different values of ry(0): 0.0, 0.5, 1.0, 2.0, 5.0, 10.0; and for 
4 different B/T ratios: 0.0, 0.1, 0.3, 0.5. The average angle to the line of sight of the 8 
independent inclination bands, was (G)= 20°, 37°, 48°, 58°, 66°, 75°, 82°, 90°. 

All these models were produced both at 4400A for the B-band, and at 9000A, for the 
I-band. We also ran some simulations, although not a complete set, at 3500A (U-band). 

4.1. Extinction properties 

Since the adopted range of parameters conforms to the one used by BFK, we can 
directly compare our results with theirs. All the results about minor and major axis 
profiles^, total extinction, face-on corrections, isophote shapes and diameters agree 

2 As defined by (9)i = J i 6sm6d6, where 6^ and Of are the limiting values for the i-th 

i 

inclination band. 

3 The shape of the major axis profiles of systems with conspicuous bulges, inclinations 
<80° and ry(0)>2, resemble closely the Type II profiles defined by Freeman (1970). Anyway, 
as noted by the same author, they often show up also in galaxies with a low dust content, 
such as SO's, and therefore their explanation seems to require peculiar stellar distributions 
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remarkably well and therefore we refer to their work (for a more thorough discussion see 
Bianchi 1995). Most of differences between our results and those of BFK are mainly 
imputable to the choice of the extinction law. While BFK use the empirical data of Rieke 
& Lebofsky (1985), that gives Ab = 1.32 Ay, we conform to the MRN and Draine & Lee 
(1984) model (Sect. |2.2.2| , Fig. [I]) and obtain Ab = 1.28 Ay- Additional differences arise 



from the finite spatial resolution of our images and from the finite resolution in inclination 
to the line of sight, as explained in appendix [A.5|. However these differences are relevant 



only when intensity varies greatly with position and inclination (the edge-on case) 

We point out that the central surface brightness of disks with large optical depth still 
varies with inclination due to stars not immersed in dust, if indeed the dust scale height 
is smaller than for stars. In addition, in the B-band, the difference between edge-on and 
face-on central surface brightness is the same for models with 7y(0) > 2, corresponding 
to 0.3 magnitudes. This common behavior makes the study of central surface brightness 
unsuitable to the determination of the optical depth of the dust disk in spiral galaxies. As 
for the influence of scattering, the main (and obvious) result is a substantial reduction of 
the actual extinction for a given ry(0). As illustrated in Fig. || such reduction is higher for 
greater ry(0) and lower inclination, with a maximum of ~ 0.3 mag for B-band models. For 
edge-on models, the difference is less than 0.1 mag, so that neglecting the effect of scattering 
when fitting the brightness profiles of highly inclined objects can be rather justified (Ohta 
& Kodaira 1995). 

EDITOR: PLACE FIGURE | HERE. 



in addition to extinction effects. 
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4-1.1. Local extinction 

The differences between models including scattering and models with absorption only 
can also be appreciated by analyzing the local extinction. In Fig. |] we plot the central 
extinction (r = 0) for models in the B band (Ab) with i = 20° (our inclination nearest to 
face-on) as a function of the central optical depth ry(0). 

EDITOR: PLACE FIGURE | HERE. 

The straight line corresponds to the Holmberg screen model, while dotted ones to sandwich 
models (Disney et al. 1989) with different £, the dust-to-stellar thickness ratio. The left 
panel is for models with absorption only, while the right panel includes the scattering: 
screen and sandwiches are the same in both panels, and only include absorption. 

Models with absorption only show that the behavior of the central extinction in the 
pure disk models (BT00) can be described locally by a sandwich with £ = 0.6, while models 
including a bulge require inverse sandwiches, that is with dust extending higher than stars 
from the equatorial plane (£ = 1.1). This is because the bulk of the emission at the center, 
that is the bulge core, is located inside the dust layer. 

In models with scattering (right panel) we note that the central extinction is reduced 
by several tenths of magnitude since the radiation has a good probability to be diffused 
by the dust, rather than absorbed. This is because, with our assumptions about the dust 
composition and size distribution, the average albedo amounts to ~ 0.5, both in the B 
and I-band. Moreover the behavior of scattering models is very different from that of 
sandwiches. 

In Fig. |5| we show the extinction at r = on the major axis, as a function of the local 
vertical optical depth, for the same models as in Figure £|. Since the bulge contribution 
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is small at this distance, all the models present nearly the same behavior, principally 
regulated by the disk. Also, models with scattering show again that sandwiches are a poor 
approximation. It is interesting to note that for ry(0) < 0.5 scattering plays such a role 
that extinction is basically absent. 

EDITOR: PLACE FIGURE | HERE. 

4-1.2. Images and color maps 

Intensity images and B-I color maps bring additional information. In Fig. |6] we show 
images in the B-band and color maps for a bulge, a disk and a BT05 model, with inclination 
82° and ry(0)= 10. Since the extinction for bulge and disk is due to the same dust disk, the 
BT05 model is simply a superposition of the first two images. We have used the same gray 
scale for the intensity images, while color maps have different scales to evidence the details. 

EDITOR: PLACE FIGURE § HERE. 

Extinction is more pronunced in the bulge because stars and dust are not homogeneously 
mixed and the dust acts almost like a screen. However, extinction also modifies the shape 
of the disk which becomes asymmetric with respect to the major axis. Color maps are 
positive everywhere, also in pure disks (intrinsic B-I=0), indicating a general reddening. 
The red color of the outer "halo" in the BT05 maps is due to the intrinsic (B-I=l) color of 
the bulge (the different appearance is due to different gray scales). 

The expected effect of scattering on colors is to render bluer the zones of an image 
where a substantial fraction of emerging radiation is contributed by forward scattering. 
This is due to the more asymmetric phase function at shorter A (Sect. |2.2.2| ). This effect 
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should be more evident in highly inclined objects and in the regions of the (dust) disk 
situated between the observer and the bright bulge core (Elvius 1956). 

Analyzing the pure bulge color maps we do note a less red lane near the center, 
coincident with the dust disk plane, but still with B-I> 1. This "blue" lane cannot be 
accounted for by scattering effects, since it is present also in models with absorption only. 
It is due, instead, to the heavy extinction of the central parts. Where this happens, the 
observed emission is mostly due to stars located out of the dust layer, and, therefore, the 
reddening is small. The reddest parts are in fact those where the local opacity (along the 
line of sight) is high but less than about 15. We note that when scattering is included, 
due to the lower effective opacity for a given Ty(0) this effect is reduced and the central 
equatorial lane actually becames redder. The disk presents this lane as well, but the 
contrast is not so strong; as an example, we show in Fig. [7| the results along the major axis 
of a pure disk with ry(0)=5 at various inclinations: the right panel is for a model with 
scattering, while the left is for pure absorption. The study of the color profile has been 
suggested by BFK as a possible tool to infer the (local) optical depth: color gradients are 
located in regions where the optical depth along line of sight is ~ 1. Our models agree with 
this result. Besides, the shape of color profiles, and then the position of gradients, is about 
the same with or without scattering; this means that pure absorption models can provide 
sufficiently good r diagnostic. 

EDITOR: PLACE FIGURE HERE. 

We have also performed simulations in U-band. Scattering effects in this wavelength 
are similar to those in the B-band and, therefore, U-B color maps are not appreciably 
different from those obtained with simpler absorption models. 



-20- 



4.2. Polarization analysis 

We now examine the polarization properties of our models. The range of parameters 
explored is the same as for the surface brightness maps. 

While the polarization introduced by a single scattering is usually conspicuous, it will 
be seen that the predicted degree of linear polarization is low (a few percent). This is due 
to three different causes: 

1. Photons emerging from the same position in a map can be either primitive, unpolarized 
photons, emitted by stars and escaped without having been intercepted by dust, 

or scattered, polarized photons. The first cause of low polarization is, therefore, 
the dilution by primitive radiation, which is higher in the regions of higher surface 
brightness, such as the central parts of both bulges and disks. 

2. If we restrict ourselves to the contribution to a certain pixel of scattered photons 
only, the resultant polarization will be maximum if the scattered photons have similar 
histories. So, the requirement for a high degree of polarization is threefold: first, the 
main contribution must come from single scattering events; second, photons must 
be scattered within a thin region along the line of sight; third, the scattering region 
must be illuminated from a narrow solid angle. The first and second requirements 
are usually met in the parameter range of our models. The third one requires a 
compact stellar distribution with respect to the dust disk and is only met by the bulge 
component. 

3. The polarization will obviously be highest if the scattering grains are identical. As 
mentioned previously our models make use of a distribution of grains in size and 
composition, and to each size and composition our code assigns, for a given scattering 
angle, a different Mueller matrix. As a result, photons emitted by the same star and 



-21 - 



scattered at the same point in the same direction may have quite different polarization 
properties, depending on the characteristics of the scattering grain. 

The maps that we will examine are obtained by averaging the images for the various 
Stokes parameters within 10x10 pixel squares. This lower resolution has been necessary to 
compensate for the lower signal to noise figures of the polarization maps. In Fig. |8| are 
collected the linear polarization maps in the B-band at different inclinations for the emission 
from the bulge only, from the the disk only, and from a composite BT05 model. The disk 
of dust is the same in all cases, namely the one with highest optical depth, ry(0) = 10. 

EDITOR: PLACE FIGURE § HERE. 

The degree of polarization comes out to be higher for larger optical depths. For the 
BT05 model of Fig. |8| at 75° inclination, the largest polarization degree is about 1.5%; 
in the same conditions but with 7y(0) = 0.5 the maximum is 0.8%. The zones of highest 
polarization are always located on, or near, the major axis. Another notable feature of our 
models is that, while increasing Ty(0) the location of the maximum polarization moves 
outward; this is illustrated by the sequences in Fig. £5[ that show the behavior of the 
polarization degree along the major axis, with galactocentric distance for some of our 
models. 

EDITOR: PLACE FIGURE | HERE. 

Clearly, for the reasons discussed above, the main contribution comes from scattering 
of the radiation from the bulge. This component has a more compact distribution with an 
emissivity highly peaked at the center, which is seen by the scatterers in the disk almost 
as a point source. The case of a pure bulge best illustrates the dilution of polarization 
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by direct radiation; the maximum degree of polarization does not vary appreciably with 
inclination. 

Pure disks show very little polarization when seen face-on; this fact is mainly due to 
the rather flat distribution of disk stars. The polarization degree increases with inclination, 
with a maximum around 50°-80° and then decreases again slightly, due to the long line of 
sight through the disk when approaching 90°. As a rule, in composite models the scattering 
of disk radiation will provide the polarization of the inner regions while the bulge will 
mainly polarize the external disk. Note how the strong polarization of the outer disk 
in case of a pure bulge 30%) is diluted by direct disk light and reduced to 1.0-1.5%. 
As a consequence, the morphology of the polarization pattern, when seen at favorable 
inclinations (50°-80°), is different in systems with different bulge luminosities: late spirals 
will have a more uniform elliptical pattern, and early ones a two lobe structure aligned 
along the major axis. 

As for the direction of the linear polarization vector, it is invariably aligned 
perpendicularly to the major axis. This means that, for the scattered field, the E 
component perpendicular to the plane of scattering Ej_, is higher than the parallel one Ey 
for most of the directions with a sizeable phase function. This is the situation commonly 
observed in scattering by small particles, where x = 2na/\ <C 1 , as it is the case, for 
example, in electron scattering. When dealing with dust grains with refractive indexes 
around one, this is not the only possible situation. The pattern can be inverted, in principle, 
at larger values of x, see for example fig. 25 of Van de Hulst (1957). We have performed a 
few simulations to check if, with larger grains or at smaller A, the pattern would appreciably 
change, in the sense of a polarization pattern with a definite alignment along the equator, 
when seen edge-on. Our results indicate that, for any reasonable dust composition and size 
distribution, this situation is difficult to achieve for A > 2000 A. The reason is that, at such 
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wavelenghts, the scattered field is mainly contributed by grafite (_L) grains which have a 
(real) refractive index m > 2. This requires, for a radical change in the pattern direction, 
values of x > 1.5, which means, in the visual (A = 5500A), a lower cutoff a_ > 0.13/im 
instead of the canonical a_ = 0.005/mi. The resulting polarization map for such a model is 
illustrated in Fig. [LC] (note the low polarization degree). 



EDITOR: PLACE FIGURE ID HERE. 



Another point to mention is the comparison between the polarization maps and the 
isophotal ones: two aspects are noticeable. First, althought the brightest, inner parts show 
very little polarization, measurable values (> 1%) are diffuse over most of the galaxy image, 
well inside the 25 B-mag arcsec -2 , and, therefore, accessible to observations (Fig. 11). 



EDITOR: PLACE FIGURE [TT] HERE. 

Second, for a given model and inclination, we find little spatial correlation between 
extinction and polarization, even in highly inclined models with large 7y(0); for example, 
the equatorial dust lanes of highly inclined disks are not strongly polarized. This is another 
signature that polarization is mainly contributed by single scattering events in the outer 
layers of the disk. 

A clear prediction of our models concerns the wavelength dependence of the degree 
of linear polarization, which we find to increase from the B to the I-band. This appears 
to be a general rule, indipendently of B/T ratios, Ty(0), and inclination. Such effect is 
ascribed to the reduced extinction in the I-band and the subsequent increase of the number 
of single scatterings. Most of the effects mentioned above are clearly depicted in Fig. [7|. 
The positions of the maxima coincide with regions where the local optical depth is ~ 1, 
similarly to what previously observed for color gradients. 



If the number of scattering is larger than one, Mie scattering can in principle introduce 
a certain degree of local circular polarization. In all our model the circular polarization 
degree is contained within the noise level and never shows at levels higher then 0.1%. 



5. Summary and conclusions 

To study the extinction and polarization properties of dusty spiral galaxies, we have 
performed extensive Monte Carlo simulations using realistic galaxy models that include 
both a bulge and a disk component, simulating different Hubble types. The scattering 
properties are calculated from Mie's theory for each dust grain; in addition, we follow the 
radiative transfer of the four Stokes parameters. The grains are assumed to be spherical, 
with a power-law size distribution, and consisting either of graphite or silicate. This enables 
us to produce both photometric and polarimetric maps of the galaxy for different optical 
depths, inclinations and B/T ratios; the main results derived from the analysis of such 
maps are summarized in the following. 

In general, we find a very good agreement with the results of BFK, who used a series of 
approximations described above, concerning luminosity profiles, total extinction, isophote 
shapes and diameters. 

We have shown that the overall effect of scattering (with respect to pure absorption) is 
to reduce the extinction for a given 7y(0). Such decrement in the effective As is particularly 
evident in the central parts of bright bulges, as shown in Fig. |], although the extinction 
fractional reduction is larger for late-type than for early-type spirals and increases with 
galactocentric distance. Sandwich models, often used in the literature, are shown to fail 
dramatically in reproducing the extinction - optical depth relation when scattering is 
properly considered, and they should not be used. From the analysis of the B-I color maps 
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we conclude that there is a net reddening of the galaxy due to dust even when scattering is 
included. In fact, the latter reduces on average the B-I color by about 0.1 mag, a mild but 
appreciable effect. Radial color gradients are found to exist in regions where the optical 
depth is of order unity. 

The degree of linear polarization produced by scattering is usually small (of the order 
of a few percent), but still practically measurable. The reasons for this behavior are: (i) 
there is a dilution effect by direct stellar (unpolarized) radiation; (ii) scattering region is 
usually illuminated from a large solid angle, thus photons with different histories enter the 
same beam, and (iii) the distribution of grain sizes and materials degrades the coherence of 
the polarization. Most of the polarization comes from light from the bulge which is single 
scattered by the outer layers of the disk. This is not surprising, of course, since the bulge 
is more compact and resembles a point source. The polarization degree increases with the 
optical depth and with inclination, flattening out above ~ 80°. The polarization pattern is 
found to be perpendicular to the major axis, independently of the morphology, inclination 
and optical depth; this result depends only on the optical properties of the dust grains, and 
can be reversed only by shifting to much larger radii the lower limit for the dust size. 

It is interesting to note that some galaxies do show polarization patterns similar to 
those of our models (Scarrott, Rolph, & Semple 1990; Fendt et al. 1995). Among them, 
late type galaxies (NGC 5907, NGC 891) show polarization perpendicular to the major 
axis on the entire extent of the disk; early type galaxies (NGC 4594, NGC 4565), instead, 
show such polarization in the outer regions of the disk, whereas the central parts have 
polarization vectors consistent with dichroic absorption, that is parallel to the major axis. 
When polarization vectors are perpendicular to the major axis, the observed degree of 
linear polarization is about one percent, in agreement with our models. 

Finally, there is little correlation between extinction and polarization degree and there 



26 



is a trend of increasing polarization from the B to the I band due to the reduced extinction 
in the I band and subsequent increase of single scatterings. 

Part of this study has been supported by the Programma Vigoni. We acknowledge 
several stimulating discussions with R.-J. Dettmar, F. Menard, E. Landi degl'Innocenti 
and M. Landolfi. One of us (S. B.) wishes to thank the Arcetri Chaos team for useful 
suggestions. 



With the Monte Carlo method each event in the photon life is determined by the 
probability distribution of a pseudo random variable, p(t), that describes the statistical 
weight of the variable t in the event, and by extracting a random number R e [0, 1]. The 
value t* associated to the random number R according to p(t) is given by the general Monte 
Carlo formula 



where t m i n and t max define the interval of possible values for t. In the following Ri will 
indicate a random number between and 1. 

As a standard reference frame we adopt a cartesian one, having the z axis along the 
galactic symmetry axis, the x and y axes on the equatorial plane. The y axis coincides with 
the nodal line (intersection between sky plane and equatorial plane). 



A. Description of the Monte Carlo code 




(Al) 
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A.l. Emission 



In our models photons are emitted by stars distributed in the galactic disk and bulge. 
The first step is to determine the coordinates of the star (i.e. the photon emission point) 
according to the star density distributions. 



Using the Jaffe law, and substituting equation (g) in equation (|A1|), one finds a simple 
relation between the radial distance and a random number Ri 



where n is the truncating factor of the Jaffe bulge and r& the scale radius. The _R 1//4 and 
Jaffe distributions are relatively similar, but the relation between r& and R e depends on 
which physical parameters, and parts of the galaxy, we intend to match best. To our 
purpose, we need the relation fitting the luminosity within a given f, since this is the 
physical variable which we will directly associate with the Monte Carlo random numbers 
(this will be the case for the disk radial distribution as well). In this sense, when uniformly 
sampled over the luminosity content, the best correspondence between the two distributions 
is found for r& = 1.16 R e . 

Then we determine the polar, 9, and azimuthal, <ft, angles in spherical coordinates. 
With spherical symmetry we have: 



A. 1.1. Bulge emission 



Ri 



(A2) 




(f)* = 2ttR 2 



0* = cos-^l - 2Ra), 



(A3) 



-28- 



where R2 and R3 are two indipendent random numbers. Thus, the emission coordinates are 



x 
Vo 



f* sin 8* cos ( 
f * sin 6* sin < 
r* cos 9*. 



In the more general case of axially symmetric ellipsoids, zq becomes 



(A4) 



zq = r*e cos 9* 



(A5) 



where e is the ratio between the minor and major axis. 



A. 1.2. Disk emission 



For an exponential disk (eq. ||), the cylindrical radius of emission r* can be found 
inverting numerically 



e r <t rdr 



max r 



e r <i rdr 



1 - [l + t*)e 
1 — (1 + m)e~ 



with t* = r* jr&. From the azimuthal simmetry around the z axis: 



2irR' 2 . 



(A6) 



(A7) 



We calculate z using 



where 



/ z * z 
Z{-)dz 
~Zmax %d 

%max y 

Z{-)dz 

%max 



R', 



Zd 



exp(- I z I /z d ), 
sech 2 (z/z d ) . 



For the exponential law, we have 



t* = sgn(R' 3 - I) In [sgn(^ - - 2R' 3 )(1 - e~ k ) + 1 



(A8) 



(A9) 



(A10) 
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for the sech 2 law, the relation is 

t* = tanh" 1 {2 [tanh(Jfe)] R' 3 - tanh(£;)} , (All) 

where t* = z*/zd and k is the vertical truncation factor. Thus 

xo = r* cos (p* 

Do = r* sin <$* (A12) 

Zq = Z*. 



A. 1.3. Direction 



If the radiation is emitted isotropically by stars, the longitude and colatitude of its 
direction are 

= 2^4, = cos- 1 (l-2i2 B ), (A13) 

from which the direction cosines are computed 

l = sin 6 cos (f> 

mo = sin 9 sin <p (A14) 
n = cos 9. 
We define the unit vector v = (lo, m , n ). 

The stellar radiation is largely unpolarized, so we will assume for the initial Stokes 
parameters 



Qo 

U 



V V y 








(A15) 



In the following the Stokes parameters will be calculated assuming as reference direction 
the projection of the z axis on the plane normal to the photon direction (for unpolarized 
radiation equation |[A15|| is true in all reference frames). 
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In summary, for the emission of every photon, we need to extract five random numbers, 
three for the coordinates of the point of emission and two for the direction of emission. 



A. 2. Diffusion point 

A photon emitted in (xQ,y ,Zo) along the direction v travels undisturbed until it 
collides with a dust grain. If undisturbed, the geometrical path within the dust distribution 
would be: 

X = Xq + l t 

y = y + m t 0<t<t T (A16) 
z = z + not 

where tr is the t-value corresponding to the point of exit from the dust disk. The total 
optical depth tt along this path is 

r A (0) ftr p(9)( r ,z) 
TT ~io ~^ dt (A17) 



r(f) 



where we made use of equation (El), and r = \/x 2 + y 2 . A fraction e~ TT of the photons will 
then emerge undisturbed from the system, while the remaining 1 — e~ TT will be absorbed or 
diffused in other directions. The point of diffusion is associated with a random number Rq, 
which determines the total optical depth at the diffusion point according to 

e^di = Rq, (A18) 



o 
or 

r = -ln(l-i? 6 ). (A19) 

If t is greater than tt, the photon will exit without scattering. This procedure, however, 
is quite inefficient: if tt is small, most of the photons will exit directly. To overcome this 
problem, we have used the method of forced scattering introduced by Cashwell & Everett 
(1959; see also Witt 1977). With this method the photon is split into two components: 
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- R e - (A20) 



one, with intensity I\ = exp(— t t ), keeps travelling undisturbed and exits the dust disk, the 
other, with intensity I 2 = 1 — exp(— r T ) is forced to scatter. Accordingly, equations (|A18| ) 
and (|A19|) become 



and 

r=-ln(l-i2 6 (l-e-^)), (A21) 

respectively. In the limit tt — > oo, r tends to the value for unforced scattering. The other 
Stokes parameters are scaled accordingly. Forcing the scattering has also the effect of 
reducing the noise in the resulting images, for the same total number of photons launches, 
but it is appreciably more time consuming. We performed a series of trial simulations using 
a point source embedded in spherical envelopes of optical opacity from 1 to 10. As a result 
of these tests, we adopted the following scheme in our code: the first scattering is always 
forced; the following ones are forced only if tt > 1. Among the various tested schemes, this 
one attains the best S/N ratio within a given computational time. 

We have used a threshold value T[ im for the optical depth: if a photon travels along a 
path with tt < r^ m it is allowed to escape without being attenuated; the adopted value 
is rum = 1CT 4 . Once r, the optical depth of the impact event, has been calculated, the 
integral in equation (|A17 ), but with t\ instead of tr, is inverted to find the value of t\, the 



geometrical length travelled before impact. The coordinates of the first scattering are then 

X\ = Xq + l()ti 

Vi = Vo + m ti (A22) 

z x = z + n ti. 
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A. 3. Choice of the dust grain 

We now need to know the geometrical and optical characteristics of the dust particle 
located at Xi,yi,Zi. The probability for a photon to collide with a spherical grain of radius 
a, refractive index m, and material i, is proportional to the product of the extinction cross 
section, Cf 1 *, and the grain number density n(a) oc a -3 5 (MRN model). The differential 
collision probability can be written as 

dPi(a) oc WiOT^Cfda, (A23) 

where Wi is the statistical weight of material i. Let us call Xi( a ) the function 

Xi(a)= f dPi(a) ae[a_,a + }.; (A24) 

Xi(a) can be normalized to unity dividing it by the value xtot(a+) = J2iXi( a +) : 

™ ^ A (A25) 

Xii a ) nas now values between and Xi( a +)- The interval [0, 1] can be divided in three 
subintervals for the three materials (Si, C||, C_L) 

[0,Xsi( a +)]> 

Xsi(a+),Xsi(a+) + Xc||( a +)] , ( A26 ) 
Xsi( a +) +Xc||(a+),1 • 

In this way, by means of a random number R 7 we choose both the grain material and size. 
The grain size is found inverting the integral 

f Xi{o)da + B i = R 7 , (A27) 

J a_ 



-33- 



where Bi is the value associated with the minimum value of each subinterval (B^ = 0, 
-Soil = Xsi( a +)i etc.). By knowing the wavelength and the material, we determine the 
complex refractive index, which, together with size, characterizes completely the optical 
properties of the grain. 



A. 4. Scattering direction and polarization transfer 

If v and v x are unit vectors defining the propagation direction before and after 
scattering, respectively, we define p^o as the unit vector indicating the projection of the 
positive z axis (the galactic symmetry axis specified by the unit vector k) onto the plane 
perpendicular to v at the scattering point. As mentioned before, p^o defines the reference 
direction for the Stokes parameters of the incoming photon. 



Two angles define vi with respect to v (Fig. |T2|) : the polar angle, 9, between v e vi 



which is the actual angle of scattering, and the azimuthal angle, 0, between p z0 and s the 
projection of vi on the plane normal to v . The vector s also lies in the scattering plane 
defined by v and vj. 



EDITOR: PLACE FIGURE IT2 HERE. 



After some algebra we derive 



Pzo = } (-kn , -m n , 1 - n„) (A28) 



Since the Mueller matrix elements are more easily defined for Stokes parameters 
referred to the scattering plane, we first have to perform a rotation about vo and change to 
Sq the reference direction. The Stokes parameters of the incoming photon referred to the 
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scattering plane are then 



( i> \ 

Qo 

U 



10 

cos 20 sin 20 

— sin 20 cos 20 

1 



( i N 

Qo 
U 

V V y 



(A29) 



With our assumptions, the phase function for polarized radiation is (van de Hulst 1957, pg. 
15) 



7nr 2 Q s 



Sii + S 



12' 



Qo 



7nr 2 Q 



2f)sca 



S u (9) + Sia (9) ( 9l C os 20 + H° sin 20 



where 



a; 



27ra 
A ' 



In 



In 



(A30) 



(A31) 



Su and S12 are two elements of the Mueller matrix, which are functions only of 9 (spherical 
grains), and C sca is the scattering cross section (analogous to C ext ). The angles 9 e are 
distributed according to equation (|A30|) . The phase function, <3>, is normalized, so that 



2- 



d6 sin 9 $ 



where $ is 



$(0) 



2?r 



rf^sin^ $(0) 



(A32) 



(A33) 



/o x 2 Q sca 
$ depends only on 9 so that 9 can be derived from the Monte Carlo method, indipendently 
both of and of the polarization of the incident radiation. The scattering polar angle, 9, 
can be calculated from a random number R 8 , inverting 



$(0) sin 9d9 = R 8 . 



(A34) 
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Substituting this value of 9 in equation (|A30|) we can find 0, from 
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R9, 



(A35) 



or 



1 

2^ 



+ 



'12 



2Sn(0 



— sin 20 + — ( 1 - 
fo fo v 



cos 20 



Rg 
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Using the values of 9 and obtained before, we compute the direction cosines of vj 

sin 



/1 



mi 



7l! 



1-ng 



— sin # 



1-ng 



—IqUq cos + m sin ) + Z cos 9 



m no cos + Iq sin ) + mo cos 9 



1 — 71q Sin COS + 72q cos 0. 
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The new direction of propagation allow us to calculate the updated Stokes parameters after 
scattering (but still referred to the scattering plane) by applying the Mueller matrix 



Q? 
u? 

I v? J 



( 



Sn S12 

S12 Sn 

s 33 

-s ? 



( v \ 
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Qo 

U'o 
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S34 

5 34 S33 j 

It will be noted that equation ( |A38| ) is different from the usual relation for the scattered 
Stokes vector, since it does not include the spherical factor l/k 2 r 2 . The standard formulae 
(e.g. Bohren & Huffman 1983, eq. [3.16]) are derived for the distribution of the scattered 
electromagnetic fields, within the frame of the electromagnetic theory. In our description, 
the scattered photon is not distributed over the entire solid angle but, instead, deflected 
onto a new single direction, and attenuated by a factor equal to the albedo (which in 
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Mie's approximation is independent of direction). If 1[ is the correct value of the scattered 
intensity, we have: 

i; = Ql' Q = Cul , (A39) 
and the correct value of the Stokes vector after scattering will be 



Q'i 



( T " \ 

o? 

u? 

V V i J 



(A40) 



The new Stokes parameters are defined with respect to s 1; the unit vector which lies in 
the scattering plane and in the plane perpendicular to the new direction of propagation vi 



(Fig. [13]): coherently with the convention adopted, we now have to calculate the Stokes 
parameters using, as a reference, p z \ the projection of the z axis onto the plane normal to 



EDITOR: PLACE FIGURE 13 HERE. 



The expression for p 2i is analogous to that for p zQ , provided we substitute vi, for v . 
Since Si lies in the scattering plane, it is parallel to the projection of vo onto the plane 
normal to vi (but with opposite direction), and we have 



si x (v - (v ■ Vi)Vi) = 



(A41) 



From the previous relation we have 
1 



si 



(/i(v • vi) - l , mi(v • vi) - m , ni(v • vi) - n ) 



'1 - (v • Vi) 2 

The angle 7 between p^i and Si is found using two other relations: 

ni(v ■ vi) - n 



(A42) 



cos 7 = p^i • Si 



1 - n(Jl - (v • v iy 



(A43) 
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sm7Vi = p z i x si 



=Vi. 



(A44) 



^l-nf^l-tvo-vi) 

The Stokes parameters referred to p^i will be calculated applying a rotation by an angle 
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In summary, every diffusion event requires four random numbers, one to locate the 
diffusion point, one for the optical characteristics of the grain and two to determine the new 
direction of propagation. 



A. 5. Imaging 

After N scatterings, the photon parameters satisfy the exit conditions. The final 
parameters are the point of last scattering, xn,Un,zn, the exit direction, VN=(lN,nT>N,nN), 
and the exit Stokes parameters, Ia^Qa^Ua^Vat. 

Our galactic models have two simmetries that can be used to reduce the computational 
time: (i) axial simmetry around the z axis and (ii) planar simmetry about the equatorial 
(xy) plane. If we look at the galaxy from a point at infinite distance in the xz plane, given 
the axial simmetry, we perform a rotation around the z axis and align v n parallel to the xz 
plane. After rotation, the exit direction can be identified by the new coordinates of the point 
of last scattering (x,y,z) and by the polar angle 6 only (we omit pedix N for simplicity). 
The corresponding Stokes parameters (I,Q,U,V) will be the same as I n ,Q„,U„,V„ because a 
rotation around the z axis does not modify the relative orientation of the wave electric field 
and the z axis. Because of the axial simmetry, for each photon exiting from (x,y,z) with 
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direction 9 and parameters (I,Q,U,V) we can add another photon exiting from (x,—y,z) with 
the same direction and parameters (I,Q,-U,-V). We can then exploit the planar simmetry: 
for each photon exiting from (x,y,z) with direction given by 9 and parameters (I,Q,U,V) 
we add a photon with (x,y,—z), direction n — 9 and parameters (I,Q,-U,-V). Using these 
simmetries we gain, therefore, a factor four in the number of photons. 

To produce maps of the galaxy as seen from different inclinations we select the photons, 
after their exit, according to their 9 direction. This is done by dividing the whole solid angle 
in iVb latitudinal bands, of equal solid angle. All the photons with exit direction within 
a given band will contribute to the image relative to that range of 9 values. Note that 
this introduces a finite resolution in inclination: for example with our standard (Nb—15) 
sampling, the most face-on image includes inclinations ranging from 0° to 30°, while the 
edge-on case includes inclinations between 86° and 94°. Having already exploited the 
planar symmetry, images with inclination 9 and n — 9 will be identical. Finally, the photons 
pertaining to a certain inclination band are projected onto the plane of the sky, according 
to their point of last scattering and direction. 

The final result consists of N B images for each Stokes parameter. Linear polarization 
maps are obtained calculating, for each pixel, the linear polarization, 

p _ VWTW 
i 

and the polarization angle, 



(A46) 



i\) = ^atan^ + (A47) 
where ipo is defined in Calamai, Landi degl'Innocenti, &Landi degl'Innocenti (1975). 
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Fig. 1. — The extinction curve resulting from the adopted dust model. Also shown are 
the contributions of the single components. The color excesses for the single materials are 
divided by the value of E(B-V) for the total mixture. 

Fig. 2. — Average phase functions as computed by the Monte Carlo code (solid lines) and 
HG phase functions with g corresponding to the mean value for the same dust distribution 
(dashed lines); the phase functions are shown at the central wavelength of the B and I 
bandpasses. 

Fig. 3. — Difference in B total magnitude between models including scattering, Bx(sc), and 
with absorption only, .By(abs), as a function of inclination, for different bulge/total ratios. 
In each panel the five lines refer to different total opacities: 7y(0) = 10, 5, 2, 1, 0.5, from 
top to bottom. 

Fig. 4. — Central (r = 0) extinction versus optical depth ry(0) for B models with i = 20° 
(solid lines); left panel is for models with absorption only, right panel for absorption + 
scattering models. Dotted lines refer to sandwich models with different dust-to-starthickness 
ratios (£). The straight solid line is for the screen model. 

Fig. 5. — As in Fig. [|, but for r = r^. 

Fig. 6. — B band images (left) and B-I color maps (right) for a bulge (top), a disk (center) 
and a BT05 model (bottom) as seen from an inclination of 82°. The dust disk is the same 
in the three cases, with 7y(0)= 10. 

Fig. 7. — Color profiles along the major axis for a BT00 model (pure disk) with ry(0) = 5 
at various inclinations. Left panel is for models with only absorption, the right one for 
absorption+scattering model. Lines are truncated at r = 4r d because of the poor S/N at 
larger radii. 
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Fig. 8. — B-band polarization maps for a bulge (top), a disk (middle), and a BT05 model 
(bottom). The optical depth is ry(0) = 10 and the inclinations 20°, 75°, 90° (from left to 
right). The dust disk is the same in the three cases. For each model, the linear polarization 
scale is given on the right. The dashed circle in the upper maps corresponds to the extent 
of the bulge. 

Fig. 9. — Linear polarization profiles along the major axis for BT05 models at inclinations 
20°, 75°, 90°, for B-band (left) and I-band (right). Profiles for r y (0) = 0.5, 2, 10 are plotted. 

Fig. 10. — B-band BT05 model seen from 75°, with t b (0) = 5 and a lower grain size cutoff 
a_ = 0.15/im; coordinates are measured in units of r d . 

Fig. 11. — B-band polarization maps for a BT05 model, with 7y(0) = 10, for inclinations 
i = 75° (left) and i = 82° (right). Isophotes up to 25 B-mag arcsec -2 are superimposed, 
spaced by 1 mag. The assumed face-on disk central brightness is 21.6 B-mag arcsec -2 
(Freeman 1970 ). 

Fig. 12. — Scattering geometry: v is the original direction of the photon and v x the direction 
after scattering; the other vectors are defined in the text. 

Fig. 13. — Scattering geometry in polar representation. 
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